package com.gitee.wsl.mathematics.function.random

import com.gitee.wsl.func.Supplier
import kotlin.jvm.JvmOverloads


/**
 * Implementation of a Halton sequence.
 *
 *
 * A Halton sequence is a low-discrepancy sequence generating points in the interval [0, 1] according to
 * <pre>
 * H(n) = d_0 / b + d_1 / b^2 .... d_j / b^j+1
 *
 * with
 *
 * n = d_j * b^j-1 + ... d_1 * b + d_0 * b^0
</pre> *
 * For higher dimensions, subsequent prime numbers are used as base, e.g. { 2, 3, 5 } for a Halton sequence in R^3.
 *
 *
 * Halton sequences are known to suffer from linear correlation for larger prime numbers, thus the individual digits
 * are usually scrambled. This implementation already comes with support for up to 40 dimensions with optimal weight
 * numbers from [
 * H. Chi: Scrambled quasirandom sequences and their applications](http://etd.lib.fsu.edu/theses/available/etd-07062004-140409/unrestricted/dissertation1.pdf).
 *
 *
 * The generator supports two modes:
 *
 *  * sequential generation of points: [.get]
 *  * random access to the i-th point in the sequence: [.skipTo]
 *
 *
 * @see [Halton sequence
 * @see [
 * On the Halton sequence and its scramblings](https://lirias.kuleuven.be/bitstream/123456789/131168/1/mcm2005_bartv.pdf)
 *
 * @since 3.3
](http://en.wikipedia.org/wiki/Halton_sequence) */
class HaltonSequenceGenerator @JvmOverloads constructor(
    dimension: Int,
    bases: IntArray = PRIMES,
    weights: IntArray = WEIGHTS,
) : Supplier<DoubleArray> {
    /** Space dimension.  */
    private val dimension: Int

    /**
     * Returns the index i of the next point in the Halton sequence that will be returned
     * by calling [.get].
     *
     * @return the index of the next point
     */
    /** The current index in the sequence.  */
    var nextIndex: Int
        private set

    /** The base numbers for each component.  */
    private val base: IntArray

    /** The scrambling weights for each component.  */
    private val weight: IntArray

    /**
     * Construct a new Halton sequence generator with the given base numbers and weights for each dimension.
     * The length of the bases array defines the space dimension and is required to be &gt; 0.
     *
     * @param dimension the space dimension
     * @param bases the base number for each dimension, entries should be (pairwise) prime, may not be null
     * @param weights the weights used during scrambling, may be null in which case no scrambling will be performed
     * @throws NullArgumentException if base is null
     * @throws OutOfRangeException if the space dimension is outside the range [1, len], where
     * len refers to the length of the bases array
     * @throws DimensionMismatchException if weights is non-null and the length of the input arrays differ
     */
    /**
     * Construct a new Halton sequence generator for the given space dimension.
     *
     * @param dimension the space dimension
     * @throws OutOfRangeException if the space dimension is outside the allowed range of [1, 40]
     */
    init {
        //NullArgumentException.check(bases)

        require(dimension in 1..bases.size){"OutOfRange Exception"}
        require(weights.size == bases.size){"Dimension Mismatch Exception"}

        /*
        if (dimension < 1 || dimension > bases.size) {
            throw OutOfRangeException(dimension, 1, PRIMES.size)
        }

        if (weights.size != bases.size) {
            throw DimensionMismatchException(weights.size, bases.size)
        }*/

        this.dimension = dimension
        this.base = bases
        this.weight = weights
        nextIndex = 0
    }

    /** {@inheritDoc}  */
    override fun get(): DoubleArray {
        val v = DoubleArray(dimension)
        for (i in 0 until dimension) {
            var index = nextIndex
            var f = 1.0 / base[i]

            val j = 0
            while (index > 0) {
                val digit = scramble(i, j, base[i], index % base[i])
                v[i] += f * digit
                index /= base[i] // floor( index / base )
                f /= base[i].toDouble()
            }
        }
        nextIndex++
        return v
    }

    /**
     * Performs scrambling of digit `d_j` according to the formula:
     * <pre>
     * ( weight_i * d_j ) mod base
    </pre> *
     * Implementations can override this method to do a different scrambling.
     *
     * @param i the dimension index
     * @param j the digit index
     * @param b the base for this dimension
     * @param digit the j-th digit
     * @return the scrambled digit
     */
    protected fun scramble(i: Int, j: Int, b: Int, digit: Int): Int {
        return if (weight != null) (weight[i] * digit) % b else digit
    }

    /**
     * Skip to the i-th point in the Halton sequence.
     *
     *
     * This operation can be performed in O(1).
     *
     * @param index the index in the sequence to skip to
     * @return the i-th point in the Halton sequence
     * @throws NotPositiveException if `index < 0`.
     */
    fun skipTo(index: Int): DoubleArray {
        require(index>=0){"NotPositive : $index"}
        /*
        if (index < 0) {
            throw NotPositiveException(index)
        }
        */

        nextIndex = index
        return get()
    }

    companion object {
        /** The first 40 primes.  */
        private val PRIMES = intArrayOf(
            2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53, 59, 61, 67,
            71, 73, 79, 83, 89, 97, 101, 103, 107, 109, 113, 127, 131, 137, 139,
            149, 151, 157, 163, 167, 173
        )

        /** The optimal weights used for scrambling of the first 40 dimension.  */
        private val WEIGHTS = intArrayOf(
            1, 2, 3, 3, 8, 11, 12, 14, 7, 18, 12, 13, 17, 18, 29, 14, 18, 43, 41,
            44, 40, 30, 47, 65, 71, 28, 40, 60, 79, 89, 56, 50, 52, 61, 108, 56,
            66, 63, 60, 66
        )
    }
}
